Single-qubit reaped quantum state tomography

Quantum state tomography is the experimental procedure of determining an unknown state. It is not only essential for the verification of resources and processors of quantum information but is also important in its own right with regard to the foundation of quantum mechanics. Standard methods have been elusive for large systems because of the enormous number of observables to be measured and the exponential complexity of data post-processing. Here, we propose a new scheme of quantum state tomography that requires the measurement of only three observables (acting jointly on the system and pointer) regardless of the size of the system. The system is coupled to a “pointer” of single qubit, and the wavefunction of the system is “reaped” onto the pointer upon the measurement of the system. Subsequently, standard two-state tomography on the pointer and classical post-processing are used to reconstruct the quantum state of the system. We also developed an efficient and scalable iterative maximum likelihood algorithm to estimate states from statistically incomplete data.


Results
Exact tomography. Consider a system of n particles, each of which has dimension d, such that the total dimension of the system is N := d n . Let {|x� |x = 0, . . . , N − 1} be the computational basis of the Hilbert space. Suppose that we have an ensemble of such systems, identically prepared in the unknown state |ψ� = N−1 x=0 |x� ψ x , with the "wavefunction" ψ x ∈ C , where C is the set of complex numbers. We assume that ψ 0 = 0 without a loss of generality (a physical state cannot be a null vector). Our proposed scheme is illustrated in the two equivalent quantum circuits in Fig. 1. We discuss these procedures in the following order: First, we select a qubit as the "pointer". The pointer plays a central role in the proposed scheme. Initially, we prepare the pointer in the state |+� := (|0� + |1�)/ √ 2 , where |0� and |1� are the computational basis states of the pointer such that the initial state of the system plus pointer is given by |�� = x |x� ψ x ⊗ |+�.
Next, we couple the system and pointer for a certain time, which is assumed to be sufficiently short compared to the typical time scales of the system and pointer. This interaction can be described by a unitary operator of the form 15 Û int = exp iθP ⊗ |1� �1| , where P is an observable of the system. For the sake of physical implementation in actual experiments, one can take two different but equivalent views of Û int . One can represent Û int with the phase shift on the pointer conditioned on the system observable P . To observe this more explicitly, let |p� be the eigenstate of the observable P belonging to the eigenvalue p and rewrite Û int as Û int = p |p� �p| ⊗Û p with the p-dependent phase shift Û p := |0� �0| + e ipθ |1� �1| on the pointer. This interpretation is depicted in the quantum circuit representation in Fig. 1a and is analogous to the conventional von Neumann picture of the measurement of the observable P . One important difference is that the pointer here is only of two dimensions and is insufficient to directly discriminate the N eigenvalues, p, of P . On the other hand, noting that Û int =Î ⊗ |0� �0| +V ⊗ |1� �1| with Î being the identity operator and V := e iθP , one can regard it as a pointercontrolled unitary operator V acting on the system. This picture is illustrated in the quantum circuit in Fig. 1b and is analogous to the quantum phase estimation circuit for a unitary transformation ( V in the present case) 16 . Throughout this paper, we will mainly consider the latter interpretation for convenience. After the unitary interaction, the total state becomes where V is the matrix representation of V in the computational basis, We then measure the eigenvalues of the observable X := x x |x� �x| in the system. When the measurement outcome is x, the (unnormalized) pointer state is reduced to Equation (3) reveals the key idea of the proposed scheme: the wavefunction ψ x appears in the two expansion coefficients and can be determined by the standard quantum state tomography by measuring three independent observables, that is, the Pauli operators σ x , σ y , and σ z in the pointer. One tricky point is that naive two-state tomography does not fix the overall phase, which is necessary to fix the relative phases of ψ x for different values of x. We now provide a careful tomographic reconstruction procedure [see Eq. (6)] that is not hindered by this tricky issue.
Physically, the two-step procedure for the measurement of X on the system and the subsequent quantum state tomography on the pointer is equivalent to the measurement of the eigenvalues of three observables, X ⊗σ z , (2) V xy := �x|V |y� = p �x|p� e ipθ �p|y� . Two equivalent schematics of the single-qubit reaped quantum state tomography. (a) The systempointer interaction is described by the p-dependent conditional phase shift Û p := |0� �0| + e ipθ |1� �1| on the pointer. (b) The system-pointer interaction is regarded as the pointer-controlled unitary operator V := e iθP on the system. The measurement on the system measures the eigenvalues x of the observable X := x x |x� �x| whereas the measurement on the pointer measures the eigenvalues of the Pauli operators σ x , σ y , or σ z . www.nature.com/scientificreports/ X ⊗σ x , and X ⊗σ y . For the purpose of mathematical analysis of measurement outcomes and maximum likelihood estimation process (see below), it is convenient to describe the measurements using the projective POVM elements

Scientific
where ˆ x = |x� �x| , ˆ m = |m� �m| , and the index m ∈ M := {0, 1, +, −, L, R} refers to the eigenstates |m� = |0� , |1� , |+� , |−� , |L� , |R� of the Pauli operators σ z , σ x , and σ y , respectively. The joint probabilities P x,m = � |Û † intˆ x,mÛint | � determine the ratio between the two coefficients, Owing to the normalization constraint, the N relations in Eq. (5) are not independent of each other. Instead of directly imposing the normalization constraint, one can just determine the ratio ψ x /ψ 0 . This casts the relation (5) to the following set of (N − 1) linear equations for x = 1, . . . , N − 1 . Given the experimentally determined measurement statistics P x,m , solving the linear equations yields the wavefunction ψ x (up to normalization). There are several dangerous cases in which Eq. (6) cannot provide a unique solution. Avoiding or overcoming them is addressed in "Methods". One simple example is to select the local basis |x� such that �x|p� = N −1/2 e 2πixk p /N , where k p is the index of p when the eigenvalues are arranged in an ordered sequence. The computational basis |x� and the eigenstates |p� of P are related by the quantum Fourier transform 17 . For a system consisting of qubits ( d = 2 ), another valuable example is the system operator of the form P = n j=1τ x j , where τ x j := (|0� �1| + |1� �0|) j denotes the Pauli operator acting on the jth qubit. This leads to a pointer-controlled unitary operator In this case, |x� and |p� are related to each other via the local Hadamard gates, with Maximum likelihood estimation algorithm. Above, we have shown that, as a matter of principle, the single-qubit reaped scheme can successfully reconstruct quantum states. It assumes an idealistic situation where the probability distribution P x,m corresponding to the POVM elements ˆ x,m can be inferred from measurements. It is possible only when the measurements are repeated infinitely many times, apart from other technical imperfections; finite repetitions give rise to statistical errors in the inferred probabilities P x,m . Obviously, the statistical errors become more severe as the system size n increases; recall the number 6d n of possible measurement outcomes (x, m). A popular method to overcome such an issue is to follow the maximum likelihood (ML) principle and seek the state that is most "likely" given the experimental observations rather than the actual (and impossible-to-infer) wavefunction [6][7][8][9]18 . In this section, we develop an iterative ML algorithm that can be combined with the single-qubit reaping scheme discussed above. We note controversies about the physically proper estimation of quantum states from the experimental data 11,18 , and it would be valuable to develop other statistical methods, such as Bayesian approaches, that are adaptable to the present tomography scheme.
Consider an ensemble of F systems. Let F x,m be the number of experimental observations corresponding to the POVM element ˆ x,m , such that F = x,m F x,m . The ideal situation corresponds to the limit F → ∞ , where the relative frequency F x,m /F gives the true probability P x,m . For finite size ( F < ∞ ), F x,m /F only estimates P x,m approximately. The observation statistics are governed by a multinomial distribution where (4) x,m :=

3ˆ
x ⊗ˆ m , www.nature.com/scientificreports/ is the probability of obtaining the result (x, m) on the condition that the system plus pointer is prepared in the state |�� = |ψ� ⊗ |+� .
We use the multinomial distribution L as the likelihood function. Generally, the likelihood function should depend on the specific measurement apparatus and other experimental conditions. Here, we focus on the generic effects on statistical error, putting aside specific technical issues. The ML approach maximizes (up to irrelevant terms) over all possible states |ψ� of the system with the normalization constraint. The wavefunction ψ x that maximizes the likelihood function satisfies the extremal equation (see "Methods" for details) where the matrix W is defined by and the x-dependent operator R x on the pointer by Formally, R x is reminiscent of a similar operator (denoted by R ) that appears in the iterative maximization algorithm adapted to the standard tomography scheme 9 . In our case, R x acts on the pointer and not on the system itself. In an ideal experiment where F → ∞ , the true wavefunction indeed gives the extremum solution, ψ x = ψ x , as R x =Î . In a realistic experiment with a finite-size ensemble ( F < ∞ ), in general ψ x � = ψ x , but ψ x is simply the wavefunction most likely for the given measurements data.
It should be noted that the operator R x depends functionally on the state |ψ� through the probability P x,m , and hence the extremum Eq. (13) is nonlinear. Solving such a nonlinear equation is unviable, particularly for large systems (involving a large number of variables ψ x ). Instead, we have developed an iterative algorithm 9,18-20 . First, we need to choose an initial trial wavefunction. From the pointer state |φ x � in Eq. (3) upon the measurement readout x, it follows that the probability P x,0 is directly proportional to |ψ x | 2 . This implies that |ψ (0) � ∝ x |x� √ F x,0 /F is a reasonable choice. At each iterative step k, the wavefunction |ψ (k) � is updated using the mapping where the iteration generator Ŵ := xy W xy |x� �y| is constructed from the matrix W in Eq. (14). Interestingly, the iteration procedure can be represented by the quantum circuit shown in Fig. 2, which illustrates the crucial role of the pointer from another perspective. The quantum circuit itself is not advantageous when one evaluates the iterations directly. However, as we will observe later, it clearly reveals the simple mathematical structure of the iteration generator Ŵ , which permits the scalability of the iterative algorithm.
The convergence of no iterative ML algorithm has been analytically proven 18 . However, in standard ML approaches 21,22 , numerical tests have demonstrated convergence for physically interesting states, and a diluted iterative algorithm is available when the convergence is critical 18 . Here, we demonstrate the algorithm numerically using several examples for a system of six qubits ( n = 6 and d = 2 ). The first example is the symmetric Dicke state |ψ� = ′ |000111� / √ 20, where ′ refers to the summation over all permutations of the qubits. We simulated the measurements for an ensemble of 24,000 systems ( F = 24, 000 ) all prepared in the same state |ψ� . The resulting relative frequencies, F x,m , of the measurement readouts (x, m) are shown in Fig. 3a. We then obtained the ML estimate |ψ (500) � for the measurement data ( F x,m ) through 500 iterations in accordance with www.nature.com/scientificreports/ (16). As shown in Fig. 3 (b, blue curve), the infidelity between the states from consecutive iterations was already less than 10 −5 after 150 iterations. The fidelity, �ψ (200) |ψ� 2 , with the true wavefunction is larger then 0.997.
We performed similar simulations and made the ML estimates for the simulation results for the W-state |ψ� = (|100000� + |010000� + · · · + |000001�)/ √ 6, the GHZ state |ψ� = (|000000� + |111111�)/ √ 2, and the ground state of the transverse-field Ising model in the ordered phase. Figure 3b corroborates the excellent convergence for all those cases. The fidelities between the ML estimates and the respective true wavefunctions were also as good as 0.99 or larger.

Scalability and mixed states. Each ML iteration in Eq. (16) involves the multiplication of exponentially
large matrices and vectors, and the computational cost of many iterations for the desired accuracy may still be high for large systems. This can be overcome by means of matrix product state (MPS) and matrix product operator (MPO) representations (see "Methods"). We first examine the quantum circuit shown in Fig. 2 more closely to better understand the MPO structure of the iteration generator, Ŵ . Let Ŵ tot be the extended operator acting on the system and pointer, which results in Ŵ = �+|Ŵ tot |+� when averaging over the pointer with the state |+� . Ŵ tot consists of the controlled-unitary operator Î ⊗ |0� �0| +V ⊗ |1� �1| and the conditional-unitary operator The former is an MPO with a bond dimension of 2 when the coupling observable P (and hence V ) is local [Eq. (7) is an example]. The latter is also an MPO with a finite bond dimension provided that the input state |ψ (k) � is an MPS with a finite bond dimension because an MPS only has finite correlations 23,24 ; see "Methods". Therefore, Ŵ tot , the product of three MPOs, should be an MPO with a finite bond dimension, and so is Ŵ as it corresponds to a partial trace of an MPO. Currently, the operation of an MPO on an MPS can be efficiently evaluated 23,24 . In summary, if the laboratory states are MPS, the iteration generator is represented by an MPO, and the ML iterations in Eq. (16) can be updated efficiently. Recently, a formally similar iterative algorithm (from a different tomography scheme) powered by MPO and MPS representations has been demonstrated in detail 22 .
Because only a polynomial number of parameters is required for the MPS representations, they span only a small portion of the entire Hilbert space. However, it is well known that many states relevant to quantum information processing, condensed matter physics, and other areas of physics exist in the MPS form. The ground states of the strongly correlated many-body Hamiltonians as well as the cluster states are notable examples.
Moreover, as was pointed out recently 14 , the tomographic estimation of MPS pure states is valuable even when the system is in a mixed state. That is, it allows us to determine a lower bound on the fidelity between the pure state estimate and mixed states compatible with the experimental observations, thereby certifying the purity of the laboratory state via experiments. A scalable ML method has been proposed to directly reconstruct mixed states via local measurements 21,22 , assuming that the states are close to a MPS. For their method, however, experimenters are required to measure many non-commuting observables whereas our scheme requires the measurement of only three observables X ⊗σ x , X ⊗σ y , and X ⊗σ z , regardless of the system size 13 .

Discussion
A seemingly similar idea to couple the system with an ancillary system and measure only one observable (over the entire system plus ancilla) has been previously proposed 25 ; this is the so-called ancilla-assisted quantum state tomography and has been demonstrated in recent experiments 26,27 . However, their scheme required the ancilla to be as large as or even larger than the system (one obvious advantage is that it can directly estimate the density matrix of the system). Moreover, no ML algorithm has been developed for their scheme.
The convergence of the ML iterations varies for different states. For example, it is noted in Fig. 3b that the convergence of the ML iterations is slower for the GHZ state (approximately 500 iterations are required for www.nature.com/scientificreports/ similar accuracy) than for other states. Recalling the massive and long distance entanglement in the GHZ state, this fact raises an interesting question about the relation between the convergence behavior of our ML iterations and the properties (such as multi-partite entanglement) of the state. We leave the relation as an inspiring open question for future works.

State-reconstruction equation.
Here, we derive the state-reconstruction Eq. (5). We begin with the (unnormalized) pointer state in Eq. (1) where we have defined α x := ψ x and β x := y V xy ψ y for notational simplicity. We want to express the ratio β x /α x in terms of the joint probabilities P x,m . The joint probabilities satisfy the following relationship: Using the last two relations, one can obtain This implies that the relative phase between α x and β x , which is the essential part for quantum coherence effects, can be extracted by combining the join probabilities on the left-hand side. More explicitly, we express it as and observe that which is identical to Eq. (5). The physical implication of the above relation is that the probabilities P x,0 and P x,1 in the computational basis of the pointer give the relative magnitudes of α x and β x , whereas the probabilities P x,± and P x,L/R give the relative phases between them.

Dangerous cases.
There are three dangerous cases where the wavefunction extraction scheme in Eq. (13) may not give a unique solution: (i) In the first case, P is compatible with the computational basis, {|x�} ( [X,P] = 0 ). Then, |x� are essentially eigenstates of P , and the pointer state upon the measurement of X becomes |φ x � = ψ x (|0� + |1� e iθx ).
Because ψ x is an overall factor, it cannot be extracted. (ii) In the second case, the unitary V is block diagonal (possibly after simultaneous permutations of rows and columns) in a given basis. Suppose that V =V (1) ⊕V (2) with V (1) and V (2) operating on orthogonal subspaces H (1) and H (2) , respectively, of H (1) ⊕ H (2) = H . Accordingly, any state |ψ� is decomposed into |ψ� = |ψ (1) � ⊕ |ψ (2) � . Upon the measurement of X , the pointer is cast to for |x� ∈ H (ν) ( ν = 1, 2 ). Therefore, in this case, one can assess ψ (ν) x /ψ (ν) 0 by applying the wavefunction extraction scheme (6) for each sector ν . However, it is impossible to extract the phase relations between different sectors. (iii) The third case is a special case where |ψ� happens to be an eigenstate of P (i.e., V ) belonging to a degenerate eigenvalue p. Suppose that the pointer is in the state |φ x � = ψ x (|0� + |1� e iθ p ) after the measurement of X on the system. The two-state tomography can successfully extract the relative phase factor e iθp , and hence p. If p is non-degenerate, the eigenvalue itself uniquely identifies |ψ� as its eigenstate. However, it is impossible if p is degenerate. Fortunately, this special case can be discerned experimentally because ϕ x is independent of x, and P x,0 = P x,1 for all x.
The first two cases can be avoided simply by properly choosing either the coupling operator P or the computational basis |x�. www.nature.com/scientificreports/ Iterative ML algorithm. Here, we detail the maximization of the likelihood function over the entire Hilbert space. Because of the normalization constraint, it is more convenient to maximize where is the Lagrange multiplier. Suppose that the system was initially in a definite state |y� and went through the unitary interaction Û int with the pointer. Let |φ xy � be the pointer state upon the measurement outcome x on the system. Explicitly, it can be expressed as The pointer state |φ x � resulting from the general initial state |ψ� of the system is related to |φ xy � by |φ x � = y |φ xy � ψ y .
In terms of |φ xy � , the joint probability can be expressed as For later use, it should be noted that its derivative with respect to ψ x has the form Then, the extremal equation for the maximization problem (23) is given by We define an x-dependent operator R x on the pointer by Then, the extremal equation (27) is Putting (24) into the above equation, we obtain which is identical to the matrix equation (13).
Matrix product states and operators. Consider a system of n particles, each of which has Hilbert space dimension d. We denote the computational basis state |x� for x = 0, 1, . . . , d n − 1 as |x� = |x 1 � ⊗ |x 2 � ⊗ · · · |x n � , where x j are the base d digits in x, x = x 1 + x 2 d + · · · + x n d n . One can observe that the conditional operator x |x� �x| ⊗R x [ψ (k) ] is an MPO with a finite bond dimension provided that the state |ψ (k) � is an MPS with a finite bond dimension. Because an MPS has finite correlations, the probabilities P x 1 ...x n ,m are factorized as they are statistically independent of the uncorrelated parts 23,24 ; we recall the base-d digits representation of x. This is also the case for the experimental observed frequencies F x 1 ...x n ,m . Therefore, the conditional operator is an MPO with a finite bond dimension.

Data availability
The source code that support the findings of this study are available from the author upon reasonable request.  (29) xz �φ xy |R x |φ xz � ψ z = ψ y .